Fractional Fokker-Planck dynamics: Numerical algorithm and simulations 
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Anomalous transport in a tilted periodic potential is investigated numerically within the frame- 
work of the fractional Fokker-Planck dynamics via the underlying CTRW. An efficient numerical 
algorithm is developed which is applicable for an arbitrary potential. This algorithm is then applied 
to investigate the fractional current and the corresponding nonlinear mobility in different wash- 
board potentials. Normal and fractional diffusion are compared through their time evolution of the 
probability density in state space. Moreover, we discuss the stationary probability density of the 
fractional current values. 
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I. INTRODUCTION 
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Thermal diffusion of Brownian particles under the 
action of a periodic force continues to present an ac- 
tive field of research over recent years, being relevant 
for various applications in condensed matter physics, 
chemicalphysics, nanotechnology, and molecular biology 

The stochastic motion of Brownian particles in a po- 
tential 



U(x) = V(x)-Fx, 



(1) 



where V(x) = V(x + L) is the periodic substrate poten- 
tial with period L and F is the constant bias, is quali- 
tatively well known Particles, subject to friction 
and noise, will diffuse and drift in the direction of the ap- 
plied bias. In the overdamped regime and in the absence 
of noise, the particles perform a creeping motion. If the 
tilting force F is large enough, so that the total poten- 
tial U(x) has no minima, the particles move down the 
corrugated plane. If minima do exist, the particles arrive 
there and stop. In the presence of noise, the particles 
do not stay permanently in the minima but will undergo 
noise activated escape events. The particles thus perform 
a hopping process from one well to the neighboring ones. 

It is well know that among many other applications 
d 0|, the model of a Brownian particle in a periodic 
potential can be used to describe Brownian motors and 
molecular motors [1,11 such as kinesins, dyneins, 

and myosins. However, many other systems, such as 
RNA polymerases, exonuclease and DNA polymerases, 
helicases, the motion of ribosomes along mRNA, the 
translocation of RNA or DNA through a pore, are ad- 
vantageously described as particles moving along a dis- 
ordered substrate. Depending on the statistical proper- 
ties of the potential, the long-time limit of the process 
can be quite different from that in a washboard poten- 
tial E3- K has been shown that the heterogeneity of 
the substrate potential may lead to anomalous dynam- 
ics [3 El 13 m In particular, over a range of forces 
around the stall force subdiffusion is observed E3i ie - 
the displacement grows as (5r 2 {t)) ~ t a , with < a < 1. 



Given the importance of the subdiffusion and the 
motion in periodic potentials in various applications 
H El E^ , and considering the biological systems men- 
tioned above, we address the physics of the effect of the 
combined action of a biased periodic force and a random 
substrate. Within our approach we model the subdiffu- 
sive dynamics in terms of a suitable residence time prob- 
ability density with a long tail El E^i El rather than 
through a random potential [l3Lll5ll20{ . 

The Fokker-Planck equation describing the over- 
damped Brownian motion in the potential U(x) can be 
generalized to anomalous transport. The corresponding 
result is known as the fractional Fokker-Planck equation 
[Til l2ll 12^ , being the central equation of fractional dy- 
namics, 



d_ 

at 



P(x,t) = D t 



d U'{x) 

dx T] a 



a dx 2 



P(x,t). (2) 



In our notation P(x, t) is the probability density, a prime 
stands for the derivative with respect to the space coordi- 
nate, K a denotes the anomalous diffusion coefficient with 
physical dimension [m 2 s~ Q ]. The quantity r] a denotes the 
generalized friction coefficient possessing the dimension 
[kgs a_2 ]; it is related to n a through r\ a K, a — k^I ', thus 
constituting a generalized Einstein relation. 

- l-a 

The notation oD t on the right-hand side of 
Eq. J2J stands for the integro-differential operator of the 
Riemann-Liouville fractional derivative, defined as fol- 
lows EHHH, 



5 l-a / x 1 

D t P(x,t) = 



At' 



P(x,f) 



(3) 



T(a)dtJ {t~t') l ~ a ' 

for < a < 1. The Riemann-Liouville operator © 
introduces a convolution integral with a slowly decay- 
ing power-law kernel, which is typical for memory effects 
in complex systems. Equation J5J describes subdiffusivc 
processes for < a < 1, and reduces to the ordinary 
Fokker-Planck equation when a = 1 . 

The fractional Fokker-Planck equation was originally 
introduced with the Riemann-Liouville fractional deriva- 
tive on its right-hand side 0, However, it can 



2 



sometimes be more convenient to switch to an equiva- 
lent representation that involves the Caputo fractional 
derivative. This formulation can provide genuine tech- 
nical advantages. The fractional Fokker-Planck equation 
can then be rewritten as: 



D?P(x,t) 



d U'{x) 
dx rj a 



a dx 2 



P(x,t), (4) 



with the Caputo fractional derivative on its left-hand 
side 0: i.e., 



1 



r(i 



at' 



i 



d 



(t - t') a dt 



7 P(x,t'). (5) 



In the present work we investigate the anomalous 
transport described by the fractional Fokker-Planck 
equation © , Q through the numerical simulation of the 
corresponding continuous time random walk (CTRW). In 
Sec.^we develop the algorithm for the numerical simula- 
tions. In Sec. IIIII we study the fractional current and the 
mobility for various types of tilted periodic potentials, 
and in Sec. lIVI the time evolution of the probability den- 
sity in space as well as the density of the current values. 
In doing so we emphasize similarities and analogies be- 
tween anomalous and normal diffusion. The implications 
of the differences are discussed in the Conclusions. 



'o='o 



O 



(I) r (2) 



(31 



FIG. 1: (color online) Sketch of the numerical algorithm: Af- 
ter a random waiting time r the particle jumps from the cur- 
rent position to the position x^ + Ax or x' n ' — Ax. 
The process is reiterated until P n ' > tfi na i. The numerical 
measurements are performed after constant time intervals at 
times t* m . The full-circles represent the events that are used 
for the computation of the physical quantities (see text for an 
explanation) . 



II. NUMERICAL SIMULATION OF THE 
FRACTIONAL FOKKER-PLANCK EQUATION 
THROUGH THE UNDERLYING CONTINUOUS 
TIME RANDOM WALK 

The fractional Fokker-Planck equation represents the 
continuous limit of a CTRW with the Mittag-Leffler res- 
idence time density plf 



E a {—{viT) a ) is the Mittag-Leffler function, 

-{VjT) a ] n 

r(na + 1) 



E a {-{ Vi TT) = Y, 



(6) 



(7) 



71=0 



and the quantity v^ 1 denotes the time-scaling parameter 
at site i. 

The numerical algorithm of the CTRW can be readily 
implemented for the motion in an arbitrary force field, 
as we will demonstrate below. 



A. Numerical algorithm for the continuous time 
random walk 



To study the CTRW in an one-dimensional potential 
U(x), we consider an ensemble of N particles moving 
on a lattice {xi = iAx}, with a lattice period Ax; i = 



0, ±1, ±2, . . . We emphasize that the numerical algorithm 
that we provide is valid for an arbitrary potential, i.e. not 
only for the potential £[]). The state of the n-th particle 
is defined through its current position x^ and the time 
t( n > at which it will perform the next jump to a nearest- 
neighbor site. 

The n-th particle of the ensemble starts from the ini- 
tial position x^fto) — x^. After a residence time r 
extracted from the probability density ipiir), the parti- 
cle jumps from site i to site i + 1 or i — 1 with proba- 
bility or q~ , respectively, obeying the normalization 
condition + gr = 1. Correspondingly, the space coor- 
dinate and the time are updated, 

x (n) _> x (n) 

± Ax and 

f( n ) — > + r. Reiterating this procedure, the full ran- 
dom trajectory of the random walker can be computed 
(see Fig. Pi. 

In order to perform the numerical simulation, one 
needs to evaluate the quantities qf 1 and Vi. In terms 
of the fractional transition rates fi and gi, they can be 
expressed in the following way, 



qt = fi/(fi + 9i), 
where we have chosen 



9i/{fi+9i): 



fi = (« Q /Ax 2 )exp[-/?([/ m - Ui)/2] 
9l = (« Q /Ax 2 )exp[-/3([/ l „ 1 - l7*)/2] 



(8) 
(9) 



(10a) 
(10b) 
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Here j3 = l/k^T is the inverse temperature and {/, = 
U(iAx). An appropriate discretisation step Ax has to 
satisfy the condition \(i(Ui±i — Ui)\ <C 1 26]. Further- 
more, the condition U"(x)Ax <C 2U'(x) must be fulfilled, 
in order to ensure the smoothness of the potential. In the 
limit Ax — > this so constructed, limiting CTRW is de- 
scribed by the fractional Fokker- Planck equation J2Jl, or 
equivalently through Eq. (0J 25j. 

In the case of a confining potential it is sufficient to 
compute the splitting probabilities qf 1 and the time scale 
parameters Vi only once over a finite x-region at the be- 
ginning of the simulation. In the case of a periodic or 
washboard potential, the quantities qf 1 and i>i can be 
computed only for the first period. In the latter case, 
while the total potential U (x) is not periodic, the poten- 
tial differences appearing in the fractional rates (|I0|l can 
be rewritten as 

U(x t ± Ax) - U{xi) = V(xi ± Ax) - V(x t ) =f FAx , 

and are therefore periodic functions of Xj. 

To perform the numerical measurements and compute 
the average {Y(t)} of a quantity Y(t) = Y(x(t)), we intro- 
duce a time lattice {t^ = mAf }, where m = 0, 1, . . . , M, 
and At* is a constant time interval between two consec- 
utive measurements. For the computation of the average 
(Y(t)}, there are at least two different strategies, which 
we discuss here. Both methods can be illustrated through 
Fig. ED 

The first possibility is as follows: Each trajectory 
x^ n \t) is separately evolved with time, until the final 
time £finai is reached, v- n ' > ifmai- As this n-th trajec- 
tory reaches a measurement time t* m (represented with 
dashed lines in Fig. QJ, i.e., > t* n , the quantity 

Ym^ — Y (x^ (t* m )) will be computed using the coordi- 
nates corresponding to the events marked with full-circles 
in Fig. ^ The value Y^fi will be saved in a storage 
variable 3^ um (t£j = S n ^m ■ After evolving all the N 
trajectories, the average is finally computed by normal- 
ization, {Y{t*J) = Y sum {t* m )/N. 

The second possibility is to evolve the whole ensemble 
until the times of all the trajectories (at which the 
particles will perform the next jump) exceed the fixed 
chosen measurement time t* m . We mark these events in 
Fig. ^ w ith full-circles. Then, all the corresponding po- 
sitions x*™) and times v- n > will be saved and the average 
(Y) at the fixed time will be computed. The proce- 
dure is reiterated to evolve the system until the final time 

^filial- 

Which of the two methods is to be preferred depends 
on the problem studied and the available computational 
resources. We use the method in which the whole ensem- 
ble is evolved in time, since it allows one to save the sys- 
tem configuration (and therefore to stop and also restart 
the time evolution) and compute the average quantities 
after each measurement time t* m . Furthermore, evolving 
the whole system together allows one to simulate a set 
of N particles interacting with each other. This method 
also allows for single-trajectory averages. 



B. Mittag-Leffler versus Pareto 

According to the Tauberian theorems [27]], for every 
< a < 1 the long time behavior of the system is de- 
termined solely by the tail of the residence time distri- 
bution 28]. Therefore, any other distribution with the 
same asymptotic form S a (i>iT) ~ l/r(I — a){viT) a could 
be used in place of the Mittag-Leffler distribution J7J). In 
fact, also the conditions S a (0) = 1 and S a (x — > oo) = 
must be satisfied, and the function S a {vir) has to de- 
crease monotonically with r. 

The Mittag-Leffler function E a {— £), defined by 
Eq. |0, can be numerically computed at £ < £o through 
the sum 



H 



h=0 



(-0 h 

r(l + ah) 



(11) 



while at values of £ > £o its asymptotic expansion can be 
used, 



K 



(0- 



^ r(l-afc) ' 



(12) 



with suitable values of H, K, and £o- 

A suitable choice for an approximate description is a 
Pareto probability density, defined by 



tf>i(r) = —7-P a {v%T) , 
dr 



with the survival probability 



[l + r(I-a) 1 /«^r] Q 



(13) 



(14) 



In the simulations of the CTRW we have usually em- 
ployed the Pareto distribution y = P a (i/jr). It is conve- 
nient numerically because it can be readily inverted to 



a 
c- 

a 

a- 1 



Stratonovich 
Pareto 

Mittag-Leffler 



-i 1 1 1 1- 




0.5 0.6 0.7 



a 



0.9 



FIG. 2: (color online). Rescaled mobility fj, a (F)ri a for 
F/Fcx = 1, as a function of a. After the same rescaled sim- 
ulation time t w 200, the distance from the asymptotic limit 
predicted by the Stratonovich formula (continuous line) of the 
mobility values corresponding the Pareto residence time den- 
sity (squares) increases significantly for a > 0.9, respect to 
those of the Mittag-Leffler residence time density (circles) . 
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provide a random residence time r [2£ 
_i y- 1/a -l 



T(l - a) 1 /" 



(15) 



y is a uniform random number in (0,1). Instead, the 
Mittag-Leffler distribution requires a specific algorithm 
to be inverted, e.g., we have used a look-up table with 
the values computed from Eqs. lfTT|l and (fT2")) . 

We have numerically verified the equivalence of the 
Mittag-Lefner and the Pareto distribution in the compu- 
tation of the asymptotic quantities. However, the dif- 
ference in the behavior of the Mittag-Leffler and the 
Pareto residence time distribution in the limit a — > 1 
has to be noticed: Namely, for a = 1 the Mittag- 
Leffler distribution transforms into the exponential func- 
tion, Ei(-UiT) = exp(— ViT), while the Pareto distribu- 
tion remains of a power-law type, leading to normal and 
anomalous diffusion, respectively. For this reason, when 
studying numerically fractional diffusion with a — > 1 
the Mittag-Leffler probability distribution should be used 
preferably. 

As the Tauberian theorems ensure the equivalence of 
the Mittag-Leffler and the Pareto distributions only in 
the asymptotic limit t — > oo, it is to be expected that at 
finite times t the two choices for the probability densities 
provide different results. The difference increases as the 
parameter a approaches the value a = 1. This situation 
is illustrated through the example in Fig. [3 



C. Summary of the algorithm 

Here we provide the core scheme of the time evolution 
algorithm used in the simulations and described above 
in Sec. Ill Al and Sec. Ill Bl The core of the program is 
the following one: 

For every measurement time t m — mAt*, where 
?7i = 1, ... , M, the loop over trajectories is performed: 

• For every trajectory n, where n = 1, . . . , N, the 
following procedure is performed: 

o While the next jumping time is smaller than 
the next measurement time, i^™' < mAt* , the 
following steps are reiterated: 

- From Eq. © the time scale parameter 
v.i at the current position i is computed. 
A random waiting time r is extracted 
from the residence time distribution, see 
Sec. Ill Bl and the next jumping time is 
computed, -> tW + r. 

- From Eqs. I© the probabilities qf to per- 
form the jump from site i to site i ± 1 are 
computed. A uniform random number be- 
tween and 1 is extracted to determine 
whether the particle jumps to the right or 



left and the new position of the particle is 
then computed, x^ — * x^ n ' ± Ax. 

o The coordinate x^ and the next jumping 
time t^ are stored. 

• Statistical averages at time t m = mAt* are com- 
puted using the stored coordinates {x^}. 



III. FRACTIONAL CURRENT AND 
GENERALIZED NONLINEAR MOBILITY IN 
WASHBOARD POTENTIALS 

Starting out with the fractional Fokker-Planck equa- 
tion in the form Q one can derive the expression for 
the mean particle position in an one dimensional tilted 
periodic potential £3> 25], reading, 



T(a + 1) 

where the stationary fractional current is given by, 

K a L [1 - exp(-/3FL)] 

~f L dxf* +L dy eM-P[U(x) - U{y)}) ' 



v a {F) 



(16) 



(17) 



This formula represents the anomalous counterpart of the 
current known for normal diffusion and reduces to the 
Stratonovich formula for a = 1 0, |3(j • For completeness 
and for the reader's convenience a simple derivation of 
Eq. i|17|) is provided in Appendix lAl 

Numerically, the fractional current is computed in the 
following manner, 



v a (F) = r(a + l) lim 

t-~ *OC 



(m - MO)) 
t a 



(18) 



The generalized nonlinear mobility is defined as 

H a (F)=v a (F)/F, (19) 

where F ^ 0. 

We test the validity of the generalized Stratonovich for- 
mula (|17fl obtained theoretically, through the simulation 
of the fractional CTRW in different periodic potentials. 
We start out with (i) the symmetric cosine potential, 



V\{x) = cos(2nx/L) 



(20) 



As another type (ii) we consider the symmetric double 
hump periodic potential 

V 2 (x) = [cos(2ttx/L) + cos(4?T2;/L)]/2 . (21) 

In order to explore the role of symmetries of the substrate 
potential we also consider (iii) the asymmetric (i.e. no re- 
flection symmetry holds), ratchet-like periodic potential, 
reading, 

V 3 (x) = [3sin(27rx/L) + sin(47nr/L)]/5 . (22) 



5 




FIG. 3: (color online). Dimensionless subcurrent 

v a {F)/(F CT /ri a ) and nonlinear mobility /j. a (F)ri a for the case 
of the cosine substrate potential (depicted in the inset) versus 
F/F CT . Numerical values corresponding to different tempera- 
tures T and fractional exponents a £ [0.1, 1] (symbols) fit the 
analytic predictions from Eq. Ijl7^ (continuous lines). 

The potentials in Eqs. J^DJl, (J2J, and as well as the 
thermal energy k^T, are measured in the same energy 
unit. For the sake of simplicity, the same symbol T is 
used in the following to represent the rescaled thermal 
energy. 

In the numerical simulations we have used the Pareto 
probability density for < a < 0.8 and the Mittag- 
Leffier density (jHJ for 0.8 < a < 1 (see the discussion in 
Sec. IIIB|) . For a = 1 corresponding to a normal Brown- 
ian process we have employed the exponential residence 
time probability density 

ipi(T) = cxp(-fjr) . (23) 

As a space step we used Ax = 0.001, measured in units 
of the space period L. For the ensemble average we 
have employed 10 4 trajectories, each one starting from 
the same initial condition x(to) — xq. The force is mea- 
sured in units of the critical tilt F CTl which corresponds 
to the disappearance of potential extrema. In the case of 
the asymmetric ratchet potential the positive critical tilt 
is used. 

We present in Figs. |3 01 andUOthe numerical results for 
the scaled fractional current v a (F)/(F CI /rj a ) and the cor- 



FIG. 4: (color online). Same as in Fig. IH1 f° r the double hump 

potential V2(x) = [cos(:r) + cos(2i)]/2. 



responding scaled nonlinear mobility; i.e. (J, a (F) rj a , with 
v a (F) and fi a (F) defined through Eq. JTSJI and Eq. i fKty . 
The subcurrent is measured in units of F cr /r] a , i.e. the 
subcurrent of a particle under the action of a constant 
bias F = F CI , the mobility is in units of the free mobility 
ry" 1 . Without loss of generality we have chosen F > 
for the symmetric substrate potentials (1201) and pijl . In 
the case of the ratchet-like potential (12211 also the results 
for negative values of the tilting force F are depicted. 

We have computed the fractional current and mobility 
for various values of a in the interval [0.1, 1]. Remarkably, 
they do not depend on the value of the fractional expo- 
nent a. For a given temperature T, all numerical values 
of Va(F)/(F CT /r)a) and ix a {F)rj a (depicted with symbols 
in Figs. 13 01 and[SJ), coincide with the theoretical curves 
resulting from Eq. I|17H (continuous lines). 

The regime of linear response at low temperatures is 
numerically not accessible. In this parameter regime the 
corresponding escape times governing the transport be- 
come far too large Q and particles are effectively trapped 
in the potential minima. At values of the tilting force 
F close to critical at which the minima disappear, the 
particles become capable to escape from the potential 
wells and the current is enhanced. The higher the tem- 
perature, the smaller is the tilting required to allow the 
particles to escape (compare the curves corresponding to 
different temperatures T in Figs. |3 01 andJSJ). At higher 
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FIG. 5: (color online). Same as in Fig. [3] for the ratchet 
potential Vk(x) = [3 sin (a; ) + sin(2x)]/5. Here also negative 
tilting is studied. 



values of the temperature T the linear response regime is 
numerically observable. 

For tilting forces F ^S> F cr or for T 3> 1 the dynamics 
approaches the behavior of a free CTRW that is exposed 
to a constant bias [l9ll25ll3l|. 



of a constant external bias, the initial probability packet 
both spreads and translates on the same time scale, one 
observes in fractional diffusion mainly a spreading only 
towards the direction of the bias [lU |31j . The probabil- 
ity density P(x,t) of particles diffusing anomalously in a 
washboard potential assumes this latter feature and, at 
the same time, undergoes the same space-periodic modu- 
lation observed in normal diffusion, being typical for mo- 
tion in a periodic potential. This is illustrated in Fig. |B| 
in which the anomalous probability density in a wash- 
board potential for a — 0.5 (lower row) is compared with 
the corresponding probability density for normal diffu- 
sion (upper row). The data of the example in Fig.[S]have 
been obtained for a tilted cosine potential with FJ F CT = 1 
and at T = 0.1. 



B. Reduced probability density 

The probability density P(x, i) associated with a nor- 
mal diffusion process in a washboard potential cannot 
relax towards a stationary, asymptotic density, due to 
the open-boundary nature of the system. However, the 
reduced asymptotic space probability density, 



p(x, t) = y2 P(nL + x,t), n 6 Z , 



(24) 



IV. 



PROBABILITY DENSITIES 



a periodic function by definition, does relax to an asymp- 
totic stationary density. Remarkably, in fractional diffu- 
sion, the probability density reaches the same stationary 
density as in the case of normal diffusion. The corre- 
sponding proof follows along the same lines of reasoning 
leading to the asymptotic fractional current v a (F) which 
is formally equivalent to the Stratonovich formula valid 
in normal diffusion [25j , as detailed also in the Appendix 
IaI This result is depicted in Fig. for the case of diffu- 
sion taking place in a tilted cosine potential. 

The form of the asymptotic reduced probability den- 
sity P st (a;) is given by 



The formal analogy between the normal and the 
fractional diffusion which emerges from the similar- 
ity between the Fokker-Planck and fractional Fokker- 
Planck equations for the current with the well-known 
Stratonovich formula and its generalization obtained for 
fractional diffusion 25] masks some basic physical differ- 
ences. For this reason we investigate and discuss here 
the time-dependent probability density in configuration 
space as well as the density of the current variable. 



A. Time evolution 

The time evolution of the space probability density 
P(x,t) in the case of anomalous diffusion is markedly 
different from that of normal diffusion (see Ref. pH l3l| ) . 
While for normal Brownian motion, under the influence 



P st (x) =M~ 1 exp[-l3U{x)] 



dx' exp[/3U{x')] 



(25) 

where M is a normalization factor (see Appendix Ia"|) . 

Even if the stationary probability density (depicted 
with continuous lines in Figs. and [SJ is the same, the 
relaxation to this stationary density is, however, very dis- 
tinct for normal and anomalous diffusion, respectively, as 
shown in Fig. [S] 

In the case of normal diffusion, at any time instant t, 
the density has only one maximum, which moves from 
the initial position (x — 0) toward its asymptotic posi- 
tion x = x'. At the same time it undergoes a spreading 
process towards the stationary density. As more particles 
reach the area around x = x' , the peak begins to grow, 
eventually spreading again to relax to the stationary so- 
lution Pst(x) (Fig. |S] above). 
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FIG. 6: Comparison between the time evolution of the probability density P(x,t) in the case of normal diffusion (above) and 
anomalous diffusion with a = 0.5 (below) in a tilted periodic potential at various evolution times t. This set up has been 
calculated for the following parameter set: k^T = 0.1, F/F cr = 1, and a cosine potential. Dotted lines at t = represent the 
initial conditions P(x,0) = S(x); i.e., all particles start out from the same position x — 0. 



In clear contrast, for a case with anomalous diffusion 
the initial probability density undergoes a spreading in 
the direction of the bias. While the initial maximum 
of the density remains at i = 0, a second maximum 
emerges at x « x', which continues to grow in weight 
as the density approaches the stationary shape P s t(x) 
(Fig. El below) . 



C. Velocity probability density 

For a particular trajectory realization x^(t), the cor- 
responding (sub)velocity reads: 



t>i n) :=T(a + l) [a?W(t)-4 n) l A" 



(26) 
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where n G l,...,N and x ( n) = x in ^(t ). This 
(sub)velocity is a random variable and one can study 
the corresponding probability density. One observes a 
spreading of the velocities corresponding to the broad 
spreading in space discussed above. The probability den- 
sity for the velocity is depicted in Fig. |5| for a periodic 
substrate cosine potential, for F/F CI — 1, T = 0.1, and 
a = 0.5. While the probability density for this velocity 
variable for normal diffusion (note the continuous line, 
right y-axis) possesses a Gaussian shape, in the anoma- 
lous case this probability density (see dashed line, left 
y-axis) assumes a very broad shape which falls off expo- 
nentially. Notably, however, the two densities have the 
same average, given by the Stratonovich formula, and 
indicated by the vertical dotted line in Fig. 



FIG. 7: (color online). Normalized theoretical stationary, 
reduced probability density P s t{x) for F/F ct = 1, T = 0.1, 
and q = 0.5, computed from Eq. (1 '_2 -~> t and l|A90 (continuous 
line) and corresponding numerical data (circles): — left y- 
axis. Also the underlying potential U(x) = cos(s) — Fx is 
depicted: — right y-axis. 



V. DISCUSSION AND RESUME 

In the field of anomalous transport the main atten- 
tion thus far has focused on the motion under the ac- 
tion of a constant or linear external force. With this 
work, we have investigated the continuous time random 
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FIG. 8: (color online) The time evolution for normal (above) 
and anomalous (below) diffusion of the reduced probability 
density P(x, i) within the first period x £ [0, L) defined 
according to Eq. 12 11 . In this example, the potential is 
U(x) — cos(x) — Fx, with F/F CT — 1, the temperature is 
T — 0.1, and the anomalous diffusion process corresponds to 
a — 0.5. The curve labels 1, 2, 3, and 4, correspond to in- 
creasing values of time; the solid line (theory) represents the 
stationary solution defined by Eq. 12511 . 



walk with power-law distributed residence times under 
the combined action of a space-periodic force and an ex- 
ternal bias, thereby elucidating in more detail parts of 
our previous findings in Ref. |25j. 

The numerical algorithm for the simulation of frac- 
tional Fokker-Planck dynamics has been detailed via the 
underlying CTRW. The application of this algorithm 
deserves to be commented on in greater detail. First, 
the effect of the replacement of the Mittag-Leffler by 
the Pareto distribution does not affect the anomalous 
transport properties in the asymptotic limit. However, 
given the finite time available for doing simulations a 
difference can still be present if the parameter a assumes 
values close to one, i.e. close to the limit of normal 
diffusion. Here, the use of the Mittag-Leffler distribu- 
tion, that precisely matches the fractional Fokker-Planck 
description, is used preferably. Otherwise, one must 
increase the overall time of simulations to arrive at 
convergent results. In order to study the fractional 
diffusion problem on the whole time scale, the use of the 
Mittag-Leffler probability density is thus unavoidable. 

Secondly, the weak ergodicity breaking [3^ makes it 
impossible to obtain the averaged value of anomalous cur- 
rent with a single time-average over a single particle tra- 
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FIG. 9: (color online) Anomalous (dashed line) and normal 
(continuous line) probability densities of the (sub)velocity, 
computed according to Eq. 12611 at t = 1000 in rescaled time 
units. The potential, tilt, temperature and a value are the 
same as in Fig. |H| The arrows point to the corresponding y- 
axes. Despite the very different shapes, the two probability 
densities possess the same average value as given by the (frac- 
tional) Stratonovich formula 1171 . indicated with the vertical 
dotted line. 



jectory. Here occurs a profound difference with the case 
of normal diffusion. Such a time-averaged quantity is it- 
self randomly distributed, as shown by the broad density 
in Fig. [5J In clear contrast to the situation with normal 
diffusion, for anomalous diffusion the current probability 
density is very broad and with a peak at the zero. Nev- 
ertheless, the average value of the current agrees very 
well with the theoretical Stratonovich value, as given by 
Eq. (|17fl . These results in turn are close in spirit to recent 
work by Bel and Barkai on the weak ergodicity breaking 
for a spatially confined fractional diffusion [33] ■ 

The results on the averaged current and its density are 
complemented by the those concerning the spatial prob- 
ability density of the particles ensemble as described by 
the fractional Fokker-Planck equation, both in the time- 
dependent case and for the appropriately chosen station- 
ary regime. The strongly asymmetric character of the 
spreading process of the probability densities in Fig. 
for the anomalous situation is clearly in line with the 
nonergodic character of the transport discussed above. 

Furthermore, a new intriguing challenge concerns the 
universality class of the single-particle current probabil- 
ity densities within the CTRW underlying the fractional 
Fokker-Planck description. The ensemble averaging is 
indispensable for obtaining a mean value of the anoma- 
lous current in agreement with the theoretical results, 
as shown for various potential shapes. In particular the 
results depicted in Fig. [S] imply the emergence of rectifi- 
cation properties when the potential tilt is alternating in 
time. The anomalous diffusion ratchet problem is, how- 
ever, highly nontrivial and complex because of intrinsic 
aging effects. For this reason, this objective is left for a 
future, separate study. We remark that it is presently not 
obvious whether an adiabatic driving limit exists at all. 
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We are confident that future work will help to clar- 
ify further and shed more light onto all these intriguing 
issues and problems. 
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APPENDIX A: STATIONARY (SUB)CURRENT 
AND PROBABILITY DENSITY IN A 
WASHBOARD POTENTIAL 

We start from the reduced probability density and the 
corresponding current, i.e., 



P(x,t) = ^P{nL + x,t), (Al) 

n 

J(x,t) = ^2j(nL + x,t) , neZ. (A2) 

n 

By definition these functions obeys periodic boundary 
conditions, P(x+L, t) = P(x, t) and J(x+L, t) = J(x, t). 
If P(x,t) is normalized, e.g. dx P(x, t) — 1, then 

P(x, t) preserves the same normalization in any x-interval 
{xq,Xq + L). The condition of stationarity, obtained by 
letting the Caputo derivative equal to zero in the conti- 
nuity equation, 



dJ(x,t) 
dx 



(A3) 



defines the reduced equilibrium probability density 
P s t(x). For both normal and fractional diffusion equa- 



tions, by integrating the resulting expression in x, one 
obtains v a /L = J st (x), i.e., explicitly, 

- V -j- = Ka exp [-(3U{x)\ — {exp \flU(x)] P st {x)} , 

(A4) 

where the integration constant v a /L is the one- 
dimensional flux. Multiplying both sides by exp [f3U(x)], 
integrating again between x and x+L, and using the con- 
ditions P Bt (x + L) = P Bt (x) and U(x + L) = U(x)-(3FL, 

pX + L 

-v a L- x j dx 1 exp \PU{x')] 

= Ka exp [l3U(x)} P st (x) [exp(-(3FL) - 1] . (A5) 

One can now multiply both sides by exp [— fiU (x)] and 
perform the final integration between x = to x = L. 
Using the normalization of P s t(^) in the x-interval (0, L), 

-v a L- x \ da; exp [-(3U(x)] / dx' exp \flU(x')] 

JO J x 

= n a [exp(-(3FL) - 1] . (A6) 

From here one can obtain the explicit expression for the 
current, 



K a L [1 - exp(-/3FL)] 



/ £ dx exp \-pU{xj\ f* +L dx' exp [f3U(x')} 



(A7) 



The stationary reduced probability density can be ob- 
tained by inverting Eq. (|A5|) . 



P st (x) = Af- 1 exp[-/3U(x)] / dx' exp[/3U(x')] ,(A8) 



x+L 



where the normalization constant is given by 

N = K a L [1 - cxp(-/3FL)] /v a . (A9) 
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